here::i_am("R/CountryVulnerabilityFigures.R")
library(here)
library(ggplot2)
library(ggridges)
library(dplyr)
library(readr)
library(sf)
library(rnaturalearth)
library(scales)
library(viridis)
library(patchwork)

# 1. Country vulnerability map---------------------------------
# Load data
country_data <- read_csv(here("data/tidy/country_vulnerability.csv"))|> 
  mutate(region = case_when(
    region == "Middle East & North Africa" ~ "Middle East & \nNorth Africa",
    region == 'East Asia & Pacific' ~ 'Latin America & \nCaribbean',
    region == 'Europe & Central Asia' ~ 'Latin America & \nCaribbean',
    region == 'Latin America & Caribbean' ~ 'Latin America & \nCaribbean',
    region == 'Sub-Saharan Africa' ~ 'Sub-Saharan \nAfrica',
    .default = region
  ))

# Get world map data
world <- ne_countries(scale = "medium", returnclass = "sf")

# Merge with vulnerability data
world_vuln <- world %>%
  left_join(country_data, by = c("iso_a3" = "country"))

# Check data range
data_range <- range(world_vuln$avg_dependency_pct, na.rm = TRUE)

# Create custom breaks based on your data
max_value <- max(country_data$avg_dependency_pct, na.rm = TRUE)
breaks <- c(0, 5, 10, 15, 20, 25, 30, 35, 40)
breaks <- breaks[breaks <= max_value]

# Add the actual maximum if it's not included
if (max_value > max(breaks)) {
  breaks <- c(breaks, ceiling(max_value))
}

# Create the map
p <- ggplot(world_vuln) +
  geom_sf(aes(fill = avg_dependency_pct), color = "white", size = 0.1) +
  # scale_fill_viridis_c(
  #   name = "Dependency\n(%)",
  #   na.value = "grey90",
  #   breaks = breaks,
  #   limits = c(0, max_value), direction = -1,
  #   labels = function(x) paste0(x, "%")
  # ) +
  scale_fill_distiller(palette = "Reds", direction = 1,
                       name = "Dependency\n(%)",
                       na.value = "grey90",
                        breaks = breaks,
                      limits = c(0, max_value),
                       labels = function(x) paste0(x, "%")) +
  theme_void() +
  theme(
    legend.position = "bottom",
    legend.key.width = unit(1.5, "cm"),
    plot.title = element_text(hjust = 0.5, size = 14),
    plot.subtitle = element_text(hjust = 0.5, size = 12, color = "grey40")
  ) +
  labs(
    title = "Global South Vulnerability to Global North Degrowth",
    subtitle = "Average dependency across all primary input categories"
  ) +
  coord_sf(crs = "+proj=robin")

# Display the plot
print(p)

# Save the plot
ggsave(here("output/VulnerabilityMap.png"), p, 
       width = 12, height = 8, 
       dpi = 300, bg = "white")

# 2. Regional vulnerability chart---------------------------------
regional_data <- read_csv("data/tidy/regional_vulnerability.csv") |> 
  mutate(region = case_when(
    region == "Middle East & North Africa" ~ "Middle East & \nNorth Africa", 
    .default = region
  )) |> 
  mutate(region = case_when(
    region == "Middle East & North Africa" ~ "Middle East & \nNorth Africa",
    region == 'East Asia & Pacific' ~ 'Latin America & \nCaribbean',
    region == 'Europe & Central Asia' ~ 'Latin America & \nCaribbean',
    region == 'Latin America & Caribbean' ~ 'Latin America & \nCaribbean',
    region == 'Sub-Saharan Africa' ~ 'Sub-Saharan \nAfrica',
    .default = region
  ))

# Create clean boxplot with jitter and mean markers
p_boxplot <- country_data %>%
  filter(region != "Other" & !is.na(avg_dependency_pct)) %>%
  ggplot(aes(
    x = reorder(region, avg_dependency_pct, FUN = median), 
    y = avg_dependency_pct, 
    fill = region)
    ) +
  geom_boxplot(
    alpha = 0.7,
    outlier.alpha = 0
  ) +
  geom_jitter(
    width = 0.2,
    color="black", size=1.5, alpha=0.5,#stroke = 0.5
  ) +
  
  # Mean markers
  stat_summary(
    fun = mean,
    geom = "point",
    shape = 23,
    size = 3,
    fill = "white",
    color = "#6F6F6F"
  ) +
  
  coord_flip() +
  scale_fill_viridis_d(option = "plasma", begin = 0.1, end = 0.9) +
  theme_minimal() +
  theme(
    legend.position = "none",
    plot.title = element_text(size = 14, hjust = 0.5),
    plot.subtitle = element_text(size = 11, hjust = 0.5, color = "grey50"),
    axis.text.y = element_text(size = 10),
    axis.text.x = element_text(size = 10),
    axis.title.y = element_blank(),
    panel.grid.major.y = element_blank(),
    panel.grid.minor = element_blank()
  ) +
  labs(
    title = "Regional Vulnerability Distribution",
    subtitle = "Box shows quartiles • Dots are countries • Diamond shows mean",
    x = "Region",
    y = "Dependency (%)",
    caption = "Global South value added dependency on Global North demand (2017)"
  ) +
  scale_y_continuous(
    labels = function(x) paste0(x, "%"),
    expand = expansion(mult = c(0.02, 0.05))
  )

# Display the plot
print(p_boxplot)

ggsave(here("output/RegionalVulnerability.png"), p_boxplot, 
       width = 12, height = 8, dpi = 300, bg = "white")

## Distribution of countries-------
country_data <- read_csv("data/tidy/country_vulnerability.csv") |> 
  mutate(region = case_when(
    region == "Middle East & North Africa" ~ "Middle East & \nNorth Africa",
    region == 'East Asia & Pacific' ~ 'Latin America & \nCaribbean',
    region == 'Europe & Central Asia' ~ 'Latin America & \nCaribbean',
    region == 'Latin America & Caribbean' ~ 'Latin America & \nCaribbean',
    .default = region
  ))

country_data_all <- country_data |> 
  mutate(region="All") |> 
  rbind(country_data)

order_regions <- c(
  'All', 
  'East Asia & \nPacific', 
  'Europe & \nCentral Asia', 
  'Latin America & \nCaribbean', 
  'Middle East & \nNorth Africa', 
  'South Asia', 
  'Sub-Saharan Africa'# , 'Other'
)
country_data_all <- country_data_all |> 
  filter(region != "Other") |> 
  mutate(region=factor(
    x=region, levels = rev(order_regions), ordered = TRUE))

p_distribution <- ggplot(
  data = country_data_all, 
  mapping = aes(
    x=avg_dependency_pct, y=region, fill=region, color=region)) +
  geom_density_ridges(
    alpha=0.5, jittered_points = TRUE,
    point_alpha=0.9, point_shape=21, point_size=1
  ) +
  scale_fill_viridis_d(
    option = "H", aesthetics=c("color", "fill")
  ) +
  theme_minimal() +
  scale_x_continuous(
    labels = label_percent(scale = 1), 
    breaks = seq(0, 50, 10)
  ) +
  labs(x = "Average dependency share", title = "Distribution of avg. dependency shares") +
  theme(
    plot.title = element_text(size = 14, hjust = 0.5),
    plot.subtitle = element_text(size = 11, hjust = 0.5, color = "grey50"),
    axis.text.y = element_text(size = 10),
    axis.text.x = element_text(size = 10),
    axis.title.y = element_blank(),
    legend.position = "none"    )

## Make new combined plot------------
library(ggpubr)
combined_plot_up <- ggarrange(p_distribution, p_boxplot, ncol = 2, labels = c("A)", "B)")) 
combined_plot_full <- ggarrange(combined_plot_up, p, nrow = 2, labels = c(NA, "C)"))

ggsave(here("output/CombinedVulnerabilityFigureFull.png"), 
       combined_plot_full, width = 14, height = 10, dpi = 300, 
       bg = "white")
